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ABSTRACT 

We  consider  the  problem  of  tracking  multiple  maneuvering  tar¬ 
gets  in  clutter  using  switching  multiple  target  motion  models. 
A  novel  suboptimal  filtering  algorithm  is  developed  by  applying 
the  basic  interacting  multiple  model  (IMM)  approach,  the  joint 
probabilistic  data  association  (JPDA)  technique  and  coupled  tar¬ 
get  state  estimation  to  a  Markovian  switching  system.  In  past 
such  an  approach  has  been  considered  using  uncoupled  target 
state  estimation.  The  algorithm  is  illustrated  via  a  simulation 
example  involving  tracking  of  two  highly  maneuvering,  at  times 
closely  spaced,  targets.  In  the  presented  example,  the  proposed 
IMM/ JPDA  coupled  filter  outperforms  an  existing  IMM/ JPDA 
uncoupled  filter. 

1.  INTRODUCTION 

We  consider  the  problem  of  tracking  multiple  maneuvering  tar¬ 
gets  in  clutter.  This  class  of  problem  has  received  considerable* 
attention  in  the  literature  [3], [6], [7],  [13],  [14].  The  switching 
multiple  model  approach  has  been  found  to  be  quite  effective 
in  modeling  highly  maneuvering  targets  [1],  [3],  [5]- [8],  [12].  In 
this  approach  various  “modes”  of  target  motion  are  represented 
by  distinct  kinematic  models,  and  in  a  Bayesian  framework,  the 
target  maneuvers  are  modeled  by  switchings  among  these  mod¬ 
els  controlled  by  a  Markov  chain.  In  the  presence  of  clutter,  the 
measurements  at  the  sensors  may  not  all  have  originated  from 
the  target-of-interest.  In  this  case  one  has  to  solve  the  problem  of 
data  association.  An  effective  approach  in  a  Bayesian  framework 
is  that  of  probabilistic  data  association  (PDA)  [3],  [5]  for  a  single 
target  in  clutter  and  that  of  joint  probabilistic  data  association 
(JPDA)  [3],  [6],  [13]  for  multiple  targets  in  clutter. 

It  is  assumed  that  the  number  of  targets  is  known  (say  N) 
and  that  for  each  target,  a  track  has  been  formed  (initiated) 
and  our  objective  is  that  of  track  maintenance.  In  [15]  such  a 
problem  has  been  considered  for  a  single  target  using  multiple 
sensors,  PDA  and  switching  multiple  models.  The  optimal  solu¬ 
tion  (in  the  minimum  mean-square  error  sense)  to  target  state 
estimation  given  sensor  measurements  and  absence  of  clutter,  re¬ 
quires  exponentially  increasing  (with  time)  computational  com¬ 
plexity;  therefore,  one  has  to  resort  to  suboptimal  approxima¬ 
tions.  For  the  switching  multiple  model  approach,  the  interact¬ 
ing  multiple  model  (IMM)  algorithm  of  [8]  has  been  found  to 
offer  a  good  compromise  between  the  computational  and  storage 
requirements  and  estimation  accuracy  [16] .  In  the  presence  of 
clutter,  one  has  to  account  for  measurements  of  uncertain  origin 
(target  or  clutter?).  Here  too,  in  a  Bayesian  framework,  one  has 
to  resort  to  approximations  to  reduce  the  computational  com¬ 
plexity,  resulting  in  the  PDA  filter  [12],  [3],  [6],  [2],  [15].  In  [15] 
the  IMM  algorithm  has  been  combined  with  the  PDA  filter  in  a 
multiple  sensor  scenario  to  propose  a  combined  IMM/MSPDAF 
(interacting  multiple  model/  multisensor  probabilistic  data  asso¬ 
ciation  filter)  algorithm.  In  [3],  [13]  and  [14]  multiple  targets  in 
clutter  (but  without  using  switching  multiple  models)  have  been 
considered  using  JPDA  filter  which,  unlike  the  PDA  filter,  ac¬ 
counts  for  the  interference  from  other  targets.  Various  versions 
of  IMMJPDA  filters  for  multiple  target  tracking  using  switching 
multiple  models  may  be  found  in  [4],  Sec.  6  of|6],  [10]  and  [11], 
While  [10]  and  [11]  present  uncoupled  filters  (i.e.  assume  that 
different  target  states  are  mutually  independent  conditioned  on 
the  past  measurements),  [4]  and  [6]  present  IMMJPDA  coupled 
filters  where  the  conditional  target  state  independence  assump¬ 
tion  is  not  made.  This  assumption  is  false  when  the  targets  are 
closely  spaced  thereby  “sharing”  measurements  (Sec.  6  of  [6]). 
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It  has  been  noted  in  [9]  that  the  IMMJPDA  coupled  filter  equa¬ 
tions  of  [4]  and  [6]  are  heuristic.  [9]  presents  an  “exact”  JPDA 
coupled  filter  for  non-switching  models  using  the  framework  of  a 
linear  descriptor  system  with  stochastic  coefficients. 

In  this  paper  we  extend  the  approach  of  [10]  which  pertains  to 
uncoupled  filtering,  to  IMMJPDA  coupled  filtering.  The  only  ap¬ 
proximations  made  are  those  typical  for  IMM  approaches;  there 
are  no  other  heuristics  as  in  [4]  and  [6].  Furthermore,  we  use  the 
standard  Markovian  switching  state-space  systems;  unlike  [9],  a 
linear  descriptor  system  framework  is  not  necessary. 

In  this  paper  we  exploit  the  basic  structure  of  [15]  in  com¬ 
bination  with  the  JPDA  algorithm  of  [3]  and  [13]  to  propose 
a  novel  IMM/ JPDA  coupled  filtering  algorithm.  As  in  [15]  it 
is  assumed  that  the  sensors  are  collocated  and  (time)  synchro¬ 
nized  with  the  sampling  rate.  Track  initialization  (formation)  is 
assumed  to  have  been  made  for  each  target.  “Standard”  assump¬ 
tions  are  used  for  JPDAF  ([13],  p.  310  of  [6]):  a  measurement 
can  have  only  one  source;  among  the  possibly  several  validated 
measurements,  at  most  one  of  them  can  be  target-originated  and 
the  remaining  validated  measurements  are  assumed  to  be  due  to 
false  alarms  or  clutter,  and  are  modeled  as  independently  and 
identically  distributed  (i.i.d.)  with  uniform  spatial  distribution 
over  the  entire  validation  region  (“across  all  targets”).  Also,  as  in 
[15]  and  [14],  we  use  sequential  updating  of  the  state  estimates 
with  measurements  (i.e.  updating  the  state  estimates  sequen¬ 
tially  with  measurements  from  different  sensors). 

2.  PROBLEM  FORMULATION 

Assume  that  there  are  total  N  targets  with  the  target  set  denoted 
as  7jv  :=  {1, 2,  •  •  • ,  N}.  Assume  that  the  dynamics  of  each  tar¬ 
get  can  be  modeled  as  one  of  the  n  hypothesized  models.  The 
model  set  is  denoted  as  M.n  ■=  {1,2,  •  •  •  ,n}  and  there  are  total 
q  sensors.  For  target  r  (r  G  7}v),  the  event  that  model  i  is  in 
effect  during  the  sampling  period  (tfc_i,£fc]  will  be  denoted  by 
M£(r).  Although  all  the  targets  share  a  common  model  set,  any 
two  targets  may  be  in  different  motion  status  from  time  to  time. 

For  the  j-th  hypothesized  model  (mode),  the  state  dynamics 
and  measurements  of  target  r  (r  G  )  are  modeled  as 

*k{r)  =  $_,(  r)*fe-i(r)  +  (,>>_,(,-)  (1) 

and 

zlk(r)  =  hi’l(xk(r))  +  w3k't(r)  for  l  =  l,--,g  (2) 

where  x k  (r)  is  the  system  state  of  target  r  at  t k  and  of  dimension 
nx  (assuming  all  targets  share  a  common  state  space),  zlk (r)  is 
the  (true)  measurement  vector  (i.e.  due  to  target  r)  from  sensor 
l  at  tk  and  of  dimension  nzU  and  Gjk_x{r)  are  the  sys¬ 

tem  matrices  when  model  j  is  in  effect  over  the  sampling  period 
(£jfc-i,£fc]  for  target  r  and  hi'1  is  the  nonlinear  transformation 
of  x*;(r)  to  zlk(r)  (l  =  1,  2 ,■••,<?)  for  model  j.  A  first-order  lin¬ 
earized  version  of  (2)  is  given  by 

4(r)  =  Hl'l(r)xk(r)  +  w^l(r)  for  1  =  1,-. ,q  (3) 

where  Hk{r')  is  the  Jacobian  matrix  of  h^'1  evaluated  at  some 
value  of  the  estimate  of  state  xk(r)  (see  Sec.  3.).  The  nature  of 
the  system  state,  the  various  matrices  in  (1)  and  (3),  and  the 
measurements  is  specified  in  more  detail  in  Sec.  4..  The  process 
noise  vje_1(r)  and  the  measurement  noise  wJkl(r)  are  mutually 
uncorrelated  zero- mean  white  Gaussian  processes  with  covari¬ 
ance  matrices  Qk_1  (same  for  all  targets)  and  R^1  (same  for  all 
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targets),  respectively.  At  the  initial  time  to,  the  initial  conditions 
for  the  system  state  of  target  r  under  each  model  j  are  assumed 
to  be  Gaussian  random  variables  with  the  known  mean  5^(r) 
and  the  known  covariance  -PpM-  The  probability  of  target  r  in 
model  j  at  to,  /^(r)  =  PfM^r)},  is  also  assumed  to  be  known. 
The  switching  from  model  Mk_x{r)  to  model  M£(r)  is  governed 
by  a  finite-state  stationary  Markov  chain  (same  for  all  targets) 
with  known  transition  probabilities  pij  —  P{M£(r)\M£_l(r)}. 
Henceforth,  tk  will  be  denoted  by  k. 

In  coupled  state  estimation  the  states  of  all  targets  are  esti¬ 
mated  jointly  [6].  To  this  end  define  the  “global”  state 

xk  :=  col{®fc(l),  Ifc(2),  •••  ,xk(N)}  (4) 

and  the  corresponding  matrices/vectors 


J  :=  col  {j  1,  •  , on}  j  jm  €  Mn  is  model  j  for  target  m, 

(5) 

FfcJ  :=  block  -  diag  {F£  (1),  Ff{ 2),  •  •  •  ,  F™  (JV)}  ,  (6) 
GJk  :=  block  -  diag  { G*  (1),  G*2  (2) ,  •  •  ■  ,  GjkN  (N) }  ,  (7) 

vk  :=  col  {i#  (1),  i42(2),  •••  ,viN(N)}.  (8) 

Then  we  have  the  state  equation  for  the  N  targets  as 

xk  =  H-\xk- 1  +  Gk-ivk-i  (9) 


The  following  notations  and  definitions  are  used  regarding  the 
measurements  at  sensor  L  Note  that,  in  general,  at  any  time  fc, 
some  measurements  may  be  due  to  clutter  and  some  due  to  the 
target,  i.e.  there  can  be  more  than  a  single  measurement  at  time 
k  at  sensor  l.  The  measurement  set  (not  yet  validated)  generated 
by  sensor  l  at  time  k  is  denoted  as 

(19) 

where  m*  is  the  number  of  measurements  generated  by  sensor  l 
at  time  k.  Variable  z1^  (i  =  1,  •  ■  •  ,m,j)  is  the  ith  measurement 
within  the  set.  The  validated  set  of  measurements  of  sensor  l  at 
time  k  will  be  denoted  by  ,  containing  m/  (<  m/)  measurement 
vectors.  The  cumulative  set  of  validated  measurements  from 
sensor  l  up  to  time  k  is  denoted  as 

=  {Yi,Y},  (20) 

The  cumulative  set  of  validated  measurements  from  all  sensors 
up  to  time  k  is  denoted  as 

Z\  =  {ZkW,Zkm,-  ■  ■  ,Zk(q)}  (21) 

where  q  is  the  number  of  sensors. 

Assuming  there  are  no  unresolved  measurements  (i.e.  mea¬ 
surement  associated  with  two  or  more  targets  simultaneously), 
any  measurement  therefore  is  either  associated  with  a  target  or 
caused  by  clutter.  Our  goal  is  to  find  the  global  state  estimate 


where  E{v£v£'}  =  Qk  :=  block  -  diag  { Q jf1 ,  •  •  •  ,  Q^N  }.  Simi¬ 
larly  define  the  global  measurement  vector  at  sensor  l  as 

4:=col{Zi(l),4(2),...,4(AT)}  (10) 


xk\k  ’■=  E{xk\ Zk} 

and  the  associated  error  covariance  matrix 

Fk\k  =  E{[xk  -  Xfc|fc][x*  -  Xk\k\'\Zi} 


(22) 

(23) 


and  the  corresponding  vectors 

hJ’l(xk)  :=  col  {V'l'W  1)),  ...  ,V'"-WA0)}  ,  (11) 

™Jkl  ■=  col  {t4> (1),  1 (2),  •  •  •  , wi» (N) }  (12) 

where  E{wk’lwk'1'}  =  Rk’1  :=  block  —  diag  j  j . 

Then  the  measurement  equation  for  N  targets  at  sensor  l  (as¬ 
suming  no  clutter  and  perfect  detections)  is  given  by 

z‘k=hJ<l(xk)+wk‘l  for  1  =  1,  •••,?.  (13) 

Following  (3),  a  first-order  linearized  version  of  (13)  is  given  by 

4  =  Hk’lxk  +  wkl  for  i  =  1,  •  |9  (14) 

where 

HJk'1  :=  diag  1),  H”'\ 2),  •••  H{n\N)}  ,  (15) 

Following  (5),  define  the  global  mode 

Mk  :=  Wk  (1).  • ' ' .  M3kN{N)}.  (16) 

The  various  targets  are  assumed  to  evolve  independently  of  each 
other.  Therefore,  the  transition  probability  for  the  global  modes 
are  given  by 


where  x ^  denotes  the  transpose  of  x k.  Included  in  the  above 
formulation  is  state  estimates  of  individual  targets. 

3.  IMM/JPDA  COUPLED  FILTERING 
ALGORITHM 

We  now  modify  the  IMM/JPDA  algorithm  of  [10]  to  apply  to 
the  coupled  system  (4)-(18);  it  will  be  called  IMM/JPDACF 
(CF  stands  for  coupled  filter).  The  approach  of  [10],  in  turn, 
is  based  on  the  approaches  of  [15],  [13],  [6]  and  [3j.  As  in  [15] 
and  [10],  for  convenience,  we  confine  our  attention  to  the  case  of 
2  sensors;  however,  the  algorithm  can  be  easily  adapted  to  the 
case  of  arbitrary  q  sensors.  As  the  IMM/MSPDAF  algorithm  is 
well-explained  in  [15]  and  Sec.  4.5  of  [6],  the  JPDAF  algorithm 
is  well-explained  in  Sec.  6.2  of  [6]  and  Sec.  9.3  of  [3]  and  the 
IMM/JPDA  filter  is  given  in  detail  in  [10]  (where  all  the  under¬ 
lying  assumptions  and  approximations  may  be  found  in  further 
detail),  we  will  only  briefly  outline  the  basic  steps  in  “one  cycle” 
(i.e.  processing  needed  to  update  for  a  new  set  of  measurements) 
of  the  IMM/JPDA  coupled  filter. 

Assumed  available:  Given  the  state  estimate  x{  ...  /  := 
E{xk-i\M^_liZ^~1}^  the  associated  covariance  and 

the  conditional  mode  probability  p>k_l  =  P[M^_1  [Zj5-1]  at  time 
k—  1  for  each  global  mode  J  G  Mn  Mn  X  ♦  •  *  X  Mn< 

Step  3.1.  Interaction  —  mixing  of  the  estimate  from  the 
previous  time  (VJ  G  Mn): 
predicted  mode  probability: 


PU  :=  P{M{'{  1),  •  •• ,  |A/*Lj(l),  •  •• ,  AC,(A0} 

N 

=  IIPi‘*  (17> 

l-l 

Similarly  we  have 

N 

d  ■=  P{M’'  (!),•••,  m=U  d  (0-  (18) 


dk~  :=  P{MJk\Zk~'}  =  Y^Pud-v  (24) 

l 

mixing  probability: 

d'J  ~  P{Ml_1\MJk,Zk1~i}  =  puni-ild'.  (25) 
mixed  estimate: 

2°kJ-i I*-,  :=  E{xk-\ | Mk ,  .2,  - 1 }  =  (26) 


covariance  of  the  mixed  estimate: 


Pfc-ijfc-i  E{[xk-i—xk_l^k_l][xk-i~xk_l\k_l]f\Mk  ,^i  } 


validation  for  each  target  separately,  the  volume  of  validation 
region  for  the  whole  target  set  is  approximated  by 


^=£^(r). 


Step  3.2.  Predicted  state  and  measurements  for  Sensor 

1  (VJeMn): 

State  prediction: 

~  E{xk\MJk,Z*-'}  =  FJk_^kJ_l[k_v  (28) 

State  prediction  error  covariance: 

p£\ ifc-i  =  £(l*fc  -*fc|fc-i][**  >-z1fc-1} 

=  Fi!-iPZ-Uk-iFi!-i  +  Gt-iQk-iGk-v  (29) 

Using  (2)  and  (28),  the  global  mode-conditioned  predicted  mea¬ 
surement  for  sensor  1  is 

zifl  *=  hJ,1(x^k_t).  (30) 

Using  the  linearized  version  (14),  the  covariance  of  the 
mode-conditioned  residual  :=  zl^  —  z^’1,  (zMJ^  := 


col{«^*l\ 


»  Zk 


}),  is  given  by 


SJkA  ■■=  + 

(31) 

where  Hk'1  is  the  first  order  derivative  (Jacobian  matrix)  of 
at  Note  that  (31)  assumes  that  zj^tr^  originates 

from  the  target  r;  the  result  (31)  does  not  depend  upon  the 
actual  measurements. 

Step  3.3.  Measurement  validation  for  sensor  1  (Vj  G 
Adn):  There  is  uncertainty  regarding  the  measurements’  ori¬ 
gins.  Therefore,  we  perform  validation  for  each  target  separately. 
There  are  two  steps  to  measurement  validation.  First  perform 
measurement  validation  for  each  target  r  (r  6  T)v)  separately. 
For  target  r,  the  validation  region  is  taken  to  be  the  same  for 
all  models,  i.e.,  as  the  largest  of  them.  Let  S3k'l(r)  denote  the 
nz i  x  nz i  submatrix  of  S^'1  including  the  rows  and  columns  of 
the  latter  numbered  as  (r  —  1)71*1  +  m,  m  =  1,2,  ••  -  ,r.  That 
is,  Sl’\r)  is  based  on  the  information  relevant  to  target  r  only. 
Let  ?kT~'l(r)  denote  the  nz\  x  1  sub-column  of  zj^’1  including  the 
rows  of  the  latter  numbered  as  (r  —  l)n*i  -f  m,  m  =  1, 2,  •  •  • ,  r; 
that  is,  zj[r,1(r)  is  the  mode-conditioned  predicted  measurement 
of  target  r  for  sensor  1.  Let  (|A|  =  det(A)) 


:=arg{  max  |S£,:l(r)l}- 
(j€Mn  J 


Then  measurement  (  i  =  1,2,  •••  ,mi)  is  validated  if  and 
only  if 

-  3r,V)Wi1  wr1!*?0  -*lr,l(r)]  <  7  (33) 

where  7  is  an  appropriate  threshold.  The  volume  of  the  valida¬ 
tion  region  with  the  threshold  7  is 

ttfto  c„,,7n*l/2|Sr1(r)|1/2  (34) 

where  nz\  is  the  dimension  of  the  measurement  and  Cnzl  is  the 
volume  of  the  unit  hypersphere  of  this  dimension  (ci  =  2,  C2  = 
7r,C3  =  47r/3,  etc.).  Choice  of  7  is  discussed  in  more  detail  in 
([6],  Sec.  2.3.2)  (see  also  Sec.  4.  later).  After  performing  the 


Step  3.4.  State  estimation  with_  validated  measure¬ 
ments  from  sensor  1  (VJ  G  Mn)-  From  among  all 
the  raw  measurements  from  sensor  1  at  time  fc,  i.e.,  Zk 

define  the  set  of  validated  measure¬ 
ment  for  sensor  1  at  time  k  as 

n1:=>fc(1).»fc(2>.-.4(™,)}  (36) 

where  ni\  is  total  number  of  validated  measurement  for  sensor  1 
at  time  k .  and 

vl(i)  ■■=  (37) 

where  1  <  h  <  fa  <  *  •  •  <  I'm  j  <  mi  when  mi  ^  0.  Note  that 
all  targets  share  a  common  validated  measurement  set  Yfc  from 
sensor  1. 

We  now  consider  joint  probabilistic  data  association  across  tar¬ 
gets  following  [6]  and  [3],  but  for  global  target  state.  A  marginal 
association  event  0*r  is  said  to  be  effective  at  time  k  when  the  val¬ 
idated  measurement  is  associated  with  (i.e.  originates  from) 
target  r  (r  =  0,  •  •  ■ ,  N  where  r  =  0  means  that  the  measurement 
is  caused  by  clutter).  Assuming  that  there  are  no  unresolved 
measurements,  a  joint  association  event  ©  is  effective  when  a 
set  of  marginal  events  {0ir}  holds  true  simultaneously.  That  is, 

©  =  n2?i0<r.  where  77  is  the  index  of  the  target  to  which  mea¬ 
surement  is  associated  in  the  event  under  consideration. 
Define  the  validation  matrix  (as  in  [6]  and  [3]) 


where  u>ir  =1  if  the  measurement  i  lies  in  the  validation  gate  of 
target  r,  else  it  is  zero.  A  joint  association  event  ©  is  represented 
by  the  event  matrix 

ft(©)  =  p*r(©)]  i  =  !,••*,  mi,  r  =  0,  •  ■  • ,  AT  (39) 


r(©)  -  |  J 


1  if  &ir  C  © 
0  otherwise. 


A  feasible  association  event  is  one  where  a  measurement  can  have 
only  one  source 

N 


]T£ir(e)  =  i  Vt, 


and  where  at  most  one  measurement  can  originate  from  a  target 


<?r(Q)  :=  ^  ^q}ir(0)  <  1  for  r  =  1,  •  •  • ,  N.  (42) 


The  above  joint  events  ©  are  mutually  exclusive  and  exhaus¬ 
tive.  Following  the  definitions  in  [6]  and  [3],  define  the  binary 
measurement  association  indicator 


Ti(©)  :=  ^Pu>ir(©),  t=  mi, 


to  indicate  whether  the  validated  measurement  is  associ¬ 
ated  with  a  target  in  event  0.  Further,  the  number  of  false 
(unassociated)  measurements  in  event  ©  is 


^(0)  =  g[i-ri(e)]. 


(44) 


and 


We  will  limit  our  discussion  to  nonparametric  JPDA  [6].  One 
can  evaluate  the  likelihood  that  the  global  mode  is  J  as 

Afc’1 

=  Pin1 1®.  m£,  z*-1}P{&\m£,  z*-1} 

© 

=  5Zp[yfcl®.  ,-Z1fc_1]P{©}  (45) 

e 

where,  as  in  (9-31)  of  [3],  the  irrelevant  conditioning  terms  have 
been  omitted  in  the  last  line  of  (45)  and  the  conditioning  on 
mi  is  implicit  in  the  event  ©.  The  second  term  (apriori  joint 
association  probabilities)  in  the  last  line  of  (45)  turns  out  to  be 
(Sec.  6.2  of  [6],  Sec.  9.3  of  [3]) 

N 

F{e}  =  ^n(PD)*'  (e)(l  -  Pd)1-6^  (46) 

S  =  1 

where  Pd  is  the  detection  probability  at  sensor  1  (assumed  to 
be  the  same  for  all  targets)  and  e  >  0  is  a  “diffuse”  prior  (for 
nonparametric  modeling  of  clutter)  whose  exact  value  is  irrele¬ 
vant.  Unlike  [10],  we  do  not  assume  that  the  states  of  the  targets 
(including  the  modes)  conditioned  on  the  past  observations  are 
mutually  independent.  Then  the  first  term  in  the  last  line  of  (45) 
can  be  written  as 

=  Vf*(e)p[Vfc(©)|Affc  ,  2f-1]  (47) 


s£'\&)  ■■=  ^{*'fc,l(©)Pfc,1(G)'|M/,2{s-1} 

=  HJk'He)P^k-iHJk'\ey  +  Hi’1.  (53) 

The  probability  of  the  joint  association  event  ©  given  that  global 
mode  J  is  effective  from  time  k  —  1  through  k  is 

/Jfc'^e)  ==  p{®\m£,z*~\y£} 

=  -cp[Y£\e,M£,zk1-1]P{e\M£,z'i-1}  =  -p[Y^MJk,z^\p{e) 

(54) 

where  the  first  term  can  be  calculated  from  (47)-(53),  the  sec¬ 
ond  term  from  (46),  and  c  is  a  normalization  constant  such  that 

Y.eP{®WJk,Zkl-\Y£}  =  l. 

Using  (from  (28))  and  its  covariance  P^k_1  (from 

(29)),  one  computes  the  partial  update  and  its  covariance 
P^k  following  the  standard  PDAF  [6],  [3],  except  that  the  global 

state  is  conditioned  on  ©,  not  the  marginal  events  0jr;  details 
follow. 

Kalman  gain:  w£(&)  =  P^_  X'WtSj^©)]-1.  (55) 

Partial  update  of  the  state  estimate: 

*$(©)  ~  E{xk\e,Mj!,Z*-\Y£} 


where  Y^©)  C  Yk  is  a  subset  of  the  validated  measurements 

consisting  of  the  measurements  associated  with  the  targets 
as  specified  by  @.  The  number  of  measurements  in  Y^©)  equal 
mj  —  <£(©)  where  <£(©)  is  the  number  of  false  alarms. 

Define  a  mi  x  [mi  —  <£(©)]  matrix  £2(9)  as  a  submatrix  of 
f)(9)  obtained  by  deleting  the  first  column  and  all  null  columns 
of  Q(9).  Then  for  a  given  ©,  we  have  a  measurement  vector 
Yk(0)  of  dimension  ri(0))n^i  given  by 

Vfc(©)  =  (/n„  t  =  1,2,  ••  •  ,mi}  (48) 

where  we  stack  up  all  target-associated  validated  measurements 
in  ©  in  ascending  order  of  targets,  In  is  the  n  X  n  identity  ma¬ 
trix,  and  the  symbol  <g)  denotes  the  Kronecker  product.  Define  a 
[(rni-^(0))n2i]  x  [Nnx]  matrix  H ^(O)  as  a  submatrix  of  H^1 
(see  (15))  obtained  by  deleting  all  i-th  block  rows  (nz\  x  .)  of 
h£'1  for  which  £,(©)  =  0.  That  is,  we  have  modified  (15)  to  keep 
only  the  block  elements  associated  with  target-associated  mea¬ 
surements  in  ©.  It  then  follows  that  the  linearized  measurement 
equation  for  Yk(G)  is  given  by 

Y^G)  =  HJk’\e)xk+wJk’1.  (49) 

Conditioned  on  the  joint  association  event  ©  and  mode  J,  the 
“coupled”  innovations  is  given  by 

,,  f  K,1  (6)  -  z1'' (©)  if  <5r(0)  =  1  for  some 

^,1(©)=<  re  (1,2 

I  0  otherwise, 

(50) 

where  Sj^f©)  is  a  subvector  of  (30)  obtained  by  deleting  all  i-th 
block  rows  (nzi  x  1)  of  for  which  £*(©)  =  0.  The  conditional 
pdf  (probability  density  function)  of  the  validated  measurements 
Yk  (©)  given  their  origins  (specified  by  ©)  and  the  global  target 
mode  J,  is  given  by 

plYtmMt'Zf-1]  =  Ar(yfc1(e);^1(e))  s£\e))  (5i) 

where 

Af(x;y,P)  :=  |27rP|-1/2exp[-i(i  - {/)'P_1(a:  -  y)\  (52) 


/  +Wk(e'>uk' *(©)  if  <5.- (©)  5^  0  for  some  1  <r  <  N 

\  if  <5r(©)  =  0  Vr  e  {1,2,  ■■■.AT}. 

(56) 

**jl  :=  E{xk\M£tZ*~\Y£}  =  ^^'1(©)?^(©)  (57) 

e 

Covariance  of  x^j*  : 

=  K |fc-i  -  Y  rt\®WJk{Q)st\Q)WJk(&)' 
©:©^©0 


+  0Jk’1(&)w£  (&ywkJ(ey 

© 


Y  ^  (®)wfc  (Q)*'*  J(©) 

.  e 

_  e 

^  J  L  «  -I 

(58) 

where  ©o  denotes  ©  for  which  <5r(©)  =  0  Vr  6  {1,2,  •  •  • ,  TV}. 
Eqn.  (58)  follows  in  a  manner  similar  to  eqn.  (3.4.2-10)  in  [6], 

Step  3.5.  The  mode-conditioned  predicted  measure¬ 
ments  for  sensor  2  (VJ  €  Mn)’-  Using  (2)  and  (57),  the 
“predicted”  measurement  for  sensor  2  is  given  by 

?£2  :=/tJ'2(£^).  (59) 


Using  the  linearized  version  (3),  the  covariance  of  the  global 
mode-conditioned  residual  :=  z 2^  —  z^’2  is  given  by 


SJk’2  :=  E{uJkM1^Jk'2^’\M£,Z^-\Y£} 


_  r/A2  nJ,l  rrJfit  ,  pJ,2 

”  Hk  ^k\kHk  +  Hk 


(60) 


where  H^'2  is  the  first  order  derivative  (Jacobian  matrix)  of 
hJ’2(.)  at  ^]‘  . 


Step  3.6.  Measurement  validation  for  sensor  2:  This  is 
similar  to  Step  3.3  where  we  replace  S^1  with  Sj^’2,  z with 

zk mi  with  m2,  Vk(r)  with  V^(r),  and  Vk  with  V£.  Details 
are  similar  to  that  in  Step  3.3,  hence  are  omitted. 


Step  3.7.  Update  with  validated  measurements  for  sen¬ 
sor  2  (VJ  €  MnY  This  step  is  similar  to  Step  3.4.  Using  the  val¬ 
idated  measurements  obtained  from  Step  3.6  and  starting  from 
x^jj.  and  P^\k  >  one  comPutes  the  final  updates  x^k  and  P£ jfc, 
and  the  likelihood 

Af  :=  (61) 

The  details  are  similar  to  that  in  Step  3.4,  hence  are  omitted. 
Step  3.8.  Update  of  mode  probabilities  (V?  £  Mn,  Vr  £ 

rNy. 

4  :=  P[^J|2f]  =  Af  (62) 

where  c  is  a  normalization  constant  such  that  =  For 

individual  targets  we  have 

Skr(r):=  P{M£{r)\Zt) 

n  n  n  n 

=  £\..  ^  pik1’",'ir-uir’ir+1',",iN  (63) 

il=l  Jr — 1=1  Jr--f  1“1  JN  =  1 

with  J  =  (ii,  •  •  • ,  in  (62). 

Step  3.9.  Combination  of  the  mode-conditioned  esti¬ 
mates  (Vr  6  Tn):  The  final  global  state  estimate  update  at 
time  k  is  given  by 

(64) 

J 

and  its  covariance  is  given  by 

pk\k  =  53  {pfc|fc  +  K|Jc  -*fc|fc][S'fe|fc  -  *fc|fcl'}  /*k-  (65) 

J 

The  state  estimate  Xfcjfc(r)  for  target  r  is  the  nx-subvector  of 
Xk\k  consisting  of  elements  (r  —  l)nx  +  m,  m  =  1, 2 

Remark  1.  Compared  to  the  uncoupled  filtering  of  [10] 
where  the  equations  are  formulated  conditioned  on  marginal  as¬ 
sociation  events  0*r,  here  we  have  conditioning  on  joint  associa¬ 
tion  events  0  for  coupled  filtering.  Eqn.  (51)  does  not  decompose 
into  the  product  of  marginal  probabilities  as  in  [10]. 

4.  SIMULATION  EXAMPLE 

We  now  consider  tracking  two  highly  maneuvering  targets  in 
clutter.  We  illustrate  the  proposed  filtering  algorithm  via  this 
example. 

The  True  Trajectory:  Target  1  starts  at  location 

[21689  10840  40]  in  Cartesian  coordinates  in  meters.  The  ini- 
tal  velocity  (in  m/s)  is  [—8.3  —  399.9  0]  and  the  target  stays  at 
constant  altitude  with  a  constant  speed  of  400m/s.  Its  trajectory 
is:  a  straight  line  with  constant  velocity  between  0  and  17s,  a 
coordinated  turn  (0.15  rad/s)  with  constant  acceleration  of  60 
m/s2  between  17  and  30s,  a  straight  line  with  constant  velocity 
between  30  and  55s,  a  coordinated  turn  (0.1  rad/s)  with  constant 
acceleration  of  40  m/s2  between  55  and  70s,  and  a  straight  line 
with  constant  velocity  between  70  and  87s.  Target  2  starts  at 
location  [30000  —  3040  40]  in  Cartesian  coordinates  in  meters. 
The  initial  velocity  (in  m/s)  is  [—382  157  0]  and  the  target  stays 
at  constant  altitude  with  a  constant  speed  of  413m/s.  Its  trajec¬ 
tory  is:  a  straight  line  with  constant  velocity  between  0  and  44s, 
a  coordinated  turn  (0.075  rad/s)  with  constant  acceleration  of  30 
m/s2  between  44  and  59s,  and  a  straight  line  constant  velocity 
between  59  and  87s. 

The  Target  Motion  Models;  These  are  patterned  after 
[15].  The  motion  models  for  the  two  targets  are  identical.  In 
each  mode  the  target  dynamics  are  modeled  in  Cartesian  coor¬ 
dinates  as  Xfc(r)  =  F(r)rcfe_i  (r)  +  G(r)vfc_i(r)  where  the  state 
of  the  target  is  position,  velocity  and  acceleration  in  each  of  the 
3  Cartesian  coordinates  (x,  y  and  z).  Thus  x^(r)  is  of  dimen¬ 
sion  9  (nx=9)  for  each  target.  Three  models  are  considered  in 
the  following  discussion.  They  are  exactly  as  in  [10]  to  which 
one  is  refered  for  more  details.  The  initial  model  probabilities 
for  two  targets  are  identical:  /ij  =  0.8,  =  0.1  and  =  0.1. 


The  mode  switching  probability  matrix  fqr  two  targets  is  also 
identical  and  is  as  in  [10]. 

The  Sensors:  Two  sensors  (we  assume  collocation,  and  time 
synchronization  of  observations,  etc.)  are  used  to  obtain  the 
measurements.  The  measurements  from  sensor  l  for  model  j  are 
zk  “  hJ,i(xfc)  +  l  =  1, 2,  reflecting  range  and  azimuth  angle 
for  sensor  1  (radar),  and  azimuth  and  elevation  angles  for  sensor 
2  (infrared).  For  further  details,  see  [10].  The  measurement  noise 
w^1  for  sensor  l  is  assumed  to  be  zero-mean  white  Gaussian  with 
known  covariances  R 1  =  diag[qr,q0i]  =  diag[400m2,  49mrad2] 
with  qa i  and  qr  denoting  the  variances  for  the  radar  azimuth  and 
range  measurement  noises,  respectively,  and  R2  =  diag[t/02,qe]  = 
diag[4mrad2,  4mrad2]  with  qa 2  and  qe  denoting  the  variances  for 
the  infrared  sensor  azimuth  and  elevation  measurement  noises, 
respectively.  Both  sensors  are  assumed  to  be  located  at  the  co¬ 
ordinate  system  origin.  The  sampling  interval  was  T  =  Is  and 
it  was  assumed  that  the  probability  of  detection  Pp  —  0.997  for 
both  sensors. 

The  Clutter:  For  generating  false  measurements  in  simula- 
tions,  the  clutter  was  assumed  to  be  Poisson  distributed  with 
expected  number  of  Ai  =  20  x  10~6/m  mrad  for  sensor  1  and 
>2  =  2  x  10~4/mrad2  for  sensor  2.  These  statistics  were  used 
for  generating  the  clutter  in  all  simulations.  However,  a  non- 
parametric  clutter  model  was  used  for  implementing  all  the  al¬ 
gorithms  for  target  tracking. 

Other  Parameters:  The  gates  for  setting  up  the  validation 
regions  for  both  the  sensors  were  based  on  the  threshold  7  =  16. 
With  the  measurement  vector  of  dimension  2,  this  leads  to  a  gate 
probability  Pq  —  0.9997  (see  p.  96  of  [6]). 

Simulation  Results:  The  results  were  obtained  from  100 

Monte  Carlo  runs.  Fig.  1  shows  the  true  trajectory  of  the  two 
targets  and  the  distance  between  the  two  targets  as  a  function 
of  time.  The  two  targets  start  out  far  apart,  move  close  to  each 
other  from  30  to  45  seconds,  and  then  move  apart  again.  Fig.  2 
shows  the  results  of  the  proposed  IMM/JPDACF  based  on  100 
runs.  Fig.  3  shows  the  results  of  the  uncoupled  IMM/JPDAF  of 
[10]  based  on  100  runs.  It  is  seen  from  Fig.  3  that  there  is  a  loss 
of  track  (target  swapping  occurs  in  2  out  of  100  runs). 
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Figure  3.  Performance  of  the  IMM/JPDAF  of  [10]  based  on 
100  runs  (read  left  to  right,  top  to  bottom):  (a)  RMSE  in  po¬ 
sition.  (b)  RMSE  in  velocity,  (c)  RMSE  in  acceleration,  (d) 
CV  model  probability  P[M£(r)|Zf]  for  r  =  1,2.  (RMSE  =  root 
mean-square  error;  CV  =  constant  velocity) 


